An Efficient Method for Fully Relativistic Simulations 
of Coalescing Binary Neutron Stars 

Walter Landry 
Physics Dept., University of Utah, SLC, UT 84112 
and Center for Radiophysics and Space Research, Cornell University, Ithaca, NY 14853 

Saul A. Teukolsky 

Center for Radiophysics and Space Research, Cornell University, Ithaca, NY 14853 



O ' 

o : 
o . 

(N ■ 

■ 



m ; 
(N '• 

> ■ 
^t- : 
o . 

o ■ 
(N ■ 

o\ : 

On ■ 

cr 

i ' 



X 



The merger of two neutron stars has been proposed as 
a source of gamma-ray bursts, r-process elements, and de- 
tectable gravitational waves. Extracting information from ob- 
servations of these phenomena requires fully relativistic sim- 
ulations. Unfortunately, the only demonstrated method for 
stably evolving neutron stars requires solving elliptic equa- 
tions at each time step, adding substantially to the computa- 
tional resources required. In this paper we present a simpler, 
more efficient method. The key insight is in how we apply nu- 
merical diffusion. We perform a number of tests to validate 
the method and our implementation. We also carry out a very 
rough simulation of coalescence and extraction of the gravi- 
tational waves to show that the method is viable if realistic 
initial data are provided. 

PACS number(s): 97.80.-d, 04.25.Dm, 04.40.Dg, 97.60. Jd 



I. INTRODUCTION 

After the discovery of the first binary neutron star sys- 
tem 1913+16, it was quickly realized that the orbit was 
decaying and the two stars would collide within a Hubble 
time. Because of the enormous available gravitational 
binding energy, neutron star-neutron star mergers like 
this became a popular mechanism for 7-ray bursts (e.g., 
jij). The merger could also eject some neutron star mate- 
rial, which led some to propose that the subsequent rapid 
decompression of the nuclear density material could cre- 
ate r-process elements |^| ||. 

However, whether or not neutron star - neutron star 
mergers can explain these phenomena is still uncertain. 
What is reasonably certain is that, if general relativity 
is correct, the coalescence will emit a gravitational wave 
signal visible at cosmological distances to detectors be- 
ing built or planned today . By comparing the signals 
at coalescence to theoretical predictions, we may be able 
to learn more about the neutron stars themselves: their 
mass, radius, and internal structure. Because neutron 
stars are such dense, relativistic objects, such a compari- 
son could also lead to new understandings in the physics 
at nuclear densities as well as the first strong field test of 
general relativity. 

When the stars collide, shocks form, black holes ap- 
pear, and in general the physics is very complicated. This 



requires a fully dynamical scheme to evolve the gravita- 
tional and matter fields. There have been a number of 
Newtonian and post-Newtonian simulations (see || and 
references therein). This approach expands the gravita- 
tional field equations around the Newtonian limit, using 
the ratio of the mass to radius, M / R, as a small param- 
eter (We use units where G = c = 1). The problem 
is that the coalescence is highly relativistic. The quan- 
tity 2M/R — 0.4 for the canonical neutron star with 
M = 1.4M Q , R = 10km. It goes up to 2M/R = 1 when 
the black hole forms, a process that Newtonian and post- 
Newtonian simulations cannot even capture. 

To treat this case, there have been some attempts to 
do a fully relativistic simulation, but with limited suc- 
cess. General relativity is a very complicated theory. 
Adding in the details of a numerical implementation 
leads to a multitude of things that can go wrong. As 
a result, self-consistent hydrodynamics-l- relativity simu- 
lations have had their greatest successes in only one or 
two dimensions [ 7[-|l2|] . These calculations have to evolve 
both the hydrodynamics (the matter) and the gravita- 
tional fields (the metric). In addition, general relativity 
allows the freedom of choosing coordinates. The calcula- 
tions have all been carried out with variations of the for- 
malism developed by Arnowitt, Deser and Misner (ADM) 
p3| . These lower dimensional calculations used a few 
tricks and applied brute force to manage long evolutions. 
In three dimensions, these techniques are either unavail- 
able or impractical. 

The neutron star problem differs from evolving black 
holes alone. In one sense, black holes are easier because 
there are no difficulties associated with the hydrodynam- 
ics or the (physically) sharp matter distribution at the 
surface. On the other hand, event horizons and singular- 
ities pose significant problems for numerical implementa- 
tions in the black hole case, which have yet to be solved 
in three dimensions. 

Undaunted, Nakamura et. al. made some pioneer- 
ing calculations of coalescing neutron stars. They came 
up against the same problem that has plagued numer- 
ical relativity until recently: growing instabilities that 
quickly crash the code. Marronetti et. al. |t4j circum- 
vented this problem by using a simplified metric. Unfor- 
tunately, using the simplified metric precludes an accu- 
rate determination of the gravitational waveform. 

Font et. al. |15| have constructed a 3D code that has 
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made some short calculations [Q. They have managed 
long hydrodynamic evolutions of single neutron stars (i.e., 
they kept the gravitational fields fixed to their analytic 
values), but a long term, self-consistent, stable evolution 
remained elusive. 

Recently, Baumgarte et. al. Jl7j modified the original 
ADM equations, trying to remove pathologies from the 
equations which might lead to instabilities. The trade- 
off is that the implementation is a little more complex. 
They applied this formalism to a number of spacetimes 
p8| , including a static neutron star. In contrast to the 
work in |jq] , they evolved the gravitational fields, but 
kept the hydrodynamic variables frozen to the analytic 
solution. They were rewarded with long term stability 
for this pseudo-dynamic problem. 

Baumgarte et. al. Jl9| applied this method to a fully 
self-consistent hydrodynamic evolution of a neutron star. 
They were able to evolve the star for several dynamical 
times, but the numerical errors eventually caused the star 
to collapse upon itself. It is also not clear whether their 
choice of coordinates will work with two coalescing stars. 
The motion of two stars around the center of mass could 
drag the coordinates along. This twisting around the 
center could lead to coordinate singularities, ruining an 
otherwise physically valid simulation. 

To alleviate this twisting, Shibata |^(| adopted a more 
robust coordinate condition (technically, approximate 
minimal distortion), and got a long term, stable, self- 
consistent, 3D relativistic simulation of coalescence. Un- 
fortunately, the coordinate condition is computationally 
expensive (it requires solving elliptic equations at every 
time step), and will significantly slow down any simula- 
tion. 

This paper presents a simpler, more efficient method 
than 1 20 1 for evolving relativistic stars. We take elements 
from many different previous investigations, add in a few 
small new contributions, and combine them into a sta- 
ble code. We use the original ADM formalism to evolve 
the metric. Although numerical relativists seem to shun 
them, we choose fully harmonic coordinates because they 
may help alleviate coordinate pathologies. There have 
been some attempts to use harmonic slicing (e.g., pl[), 
but not the full harmonic gauge. We use the relatively 
new and sophisticated high resolution shock capturing 
method for relativistic hydrodynamics developed in [ p"5| . 
For the outer boundaries, we adopt the condition from 
jl7j that allows gravitational waves to propagate off the 
grid. Our new contribution is how we add in a small nu- 
merical diffusion to stabilize the simulation. The usual 
way to add it in diffuses mostly along coordinate axes 
(x, y, z). We use a new scheme that also diffuses along 
diagonals (x = y, x = —y, x = z, x = —z, y = z, 
y = —z). This small variation allows long term, accurate 
evolutions. 

We then present a number of tests: short term tests to 
validate our implementation of the equations, and long 
term tests to validate the whole approach. Then we pro- 
ceed to compute a very unrealistic simulation of coalesc- 



ing binary neutron stars to demonstrate that the method 
will work. Finally, we describe what remains to be done 
to turn this into a realistic simulation. We expect that 
reasonably accurate simulations will become available in 
the next few years, although they may require the largest 
supercomputers . 



II. METHODS 

A. Metric Evolution 

To evolve the metric, we use the standard decompo- 
sition of Arnowitt, Deser and Misner (ADM) Jlj|. The 
line element then takes the form 



ds 2 
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dx»dx v 



-a 2 dt 2 



j i3 (dx l + (3 i dt)(dx j + (3 j dt). 

(1) 



The lapse and shift functions a and /J; embody the gauge 
freedom of general relativity (i.e., the freedom to choose 
arbitrary coordinates), and so can be chosen without re- 
striction. The embedding of the t = constant slices in 
the spacetime is described geometrically by the extrinsic 
curvature K%j, which is defined by the equation 



dt 



-2aK i3 + ViPj + Vjft 



(2) 



where Vi denotes a covariant derivative with respect to 
the three dimensional metric 7^. We have written this 
equation in the form of an evolution equation for 7^ . 
Using Einstein's equations one can derive the evolution 
equation for Kij: 



dKj. 
dt 



= -V 3 a + KijV if}' 



Rij - 2K a K) 



K Kij S-ij 



(3) 



This in turn introduces the Ricci tensor R^j and the 
matter terms p and S. Rij is the three-dimensional Ricci 
tensor associated with the metric 7y . The most common 
way to write it involves first constructing the connection 
coefficients T l jk with the equation 
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■ Ilk,, 



Ijk.l), 



where commas denote partial derivatives. Then we could 
use the standard textbook formula to construct Rij . Un- 
fortunately, this involves taking derivatives of r* fe . In 
taking a derivative of a derivative numerically, there is 
no simple way to keep the error in the derivatives second 
order convergent and have the finite difference stencil in- 
clude only nearest neighbor couplings. Second order ac- 
curacy requires knowledge of points two grid points away. 
This complicates boundary conditions, because we would 
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have to apply them to a layer of points two points deep 
instead of one point deep. So, we write the Ricci tensor 
in an alternative form as 

Rij — 2^ Tiijfej TfcMi Tij,fei 

+2 (n?r mfci - r^r rofcl )] . (4) 

Now, we can take centered, second-order differences of 
the metric 7^ , for which we only need to use the nearest 
neighbors. 

p and S are projections of the four-dimensional stress 
energy tensor onto the three-dimensional spacetime 
described by 7^. They are defined in terms of the 4- 
vector normal to the three dimensional slice 

n p = (-a, 0,0,0). 

We then define three projections, p, J 1 , and Sij as 

p = Sim^T^ = 8na 2 T u , 
J 1 = -87rn M 7jT«, 

and 



B. Constraints 

The evolution equations (eqs. ^ and ^) involve only 6 
of the 10 Einstein equations. The remaining equations 
are the energy, or Hamiltonian, constraint 

R + K 2 - KijK ij = 2p, (5) 

and the three momentum constraints 

Vj (K ij - j ij K) = J\ (6) 

If we start with data that satisfies these constraints and 
evolve them with eqs. || and ||, then we are guaranteed to 
evolve to a spacetime that still satisfies the constraints. 
(Note that it doesn't matter how a and /3 Z are evolved. 
Their evolution is determined by the gauge, and eqs. |5| 
and |^ are true in any gauge.) 

This guarantee is analytic, but errors in a numerical 
evolution can lead to a (7^, Kij) pair that does not 
exactly satisfy the constraints. Accordingly, we gener- 
ally use the constraints as a check on the accuracy of 
the code. However, in order to get self-consistent initial 
data for two neutron stars, we need to explicitly enforce 
the constraints. In doing this, an ambiguity immediately 
presents itself. Eqs. [| and || do not declare what com- 
bination of the twelve components of 7^ and K^ should 
be constrained by the four constraints. A natural way, 
but by no means the only way, is to decompose 7^ and 



Kij with the York decomposition p2[ . We decompose 
the metric by defining a conformal metric 7^ = 7i J -'0 — 4 , 
where detTy = 1. The inverse of 7^ is then -y 1 - 7 = j^ip 4 . 
We then apply the energy constraint to ip. 

Decomposing the extrinsic curvature is a little more 
complicated. We decompose the extrinsic curvature into 
traced and trace-free parts 

K ij = t/T 10 + (/A) u ) + i?/r 4 7 2fe TrK, 

where 

[lX) tJ = V'X 3 + V J A* - -f-? V fe A fe , 

o 

and all covariant derivatives V are with respect to the 
conformal metric 7^ . This is not a unique decomposition, 
since and (ZA) y describe overlapping degrees of free- 
dom. In practice, we decompose K^ assuming X 1 = 0, 
solve for a new X % which makes K# solve the constraints, 
and recompose K^. 

In terms of these variables, the constraints become 

- 8V 2 r/> = -Rip - \ (trK) 2 ifi 5 

+ + (IXfY ^ + 2p^ 5 , (7) 

v 2 a i + -v i v j x j + mx j 
3 J 

= J l V 10 - VjA* + ^j 6 VHrK 7 (8) 

where R is the Ricci tensor associated with the confor- 
mal metric 7. To solve these equations, we first linearize 
them, writing ip as ipQ + 5i/j and X 1 as Aq + <5A\ We 
explicitly decouple the linearized equations by setting 
8X l = in the linearized form of eq. |t], and Sip = 
in the linearized form of eq. ||. The linearized equations 
then become 

- 8V 2 (^ + H) = -%0 + fy) - I (trK) 2 ^(4, + 5#) 
+ (l« + (lX ) ij ) 2 VV 8 (^o - 7Sip) + 2p^ Q + 55il>), (9) 

V 2 (A* + SX l ) + Jv'V^Xg + SX) + RAXl + 8X l ) 

= J'Vo - VjA iS + \ipffiHrK. (10) 

Then we solve the linearized equations using multigrid 
methods |2^|. We use red-black Gauss-Seidel smooth- 
ing, with 2 pre- and post-smoothings for the linearized 
ip equation, and 5 pre- and post-smoothings for the lin- 
earized X 1 equation. On the outer boundaries, we as- 
sume that ip = ^original and X 1 = A* rlglna] = (i.e., 
we freeze the boundaries). When working on parallel 
machines, we smooth the solution on each machine indi- 
vidually, and then synchronize the edges. Then we add 
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Sip to ipo and 8X l to Xq, recompute the terms involving 
ip Q and Xq in eqs. § and and iterate until we reach a 
tolerance of 10 -10 in the norm of eqs. ^ and [l0| 

The explicit decoupling helps the elliptic solver con- 
verge. Decoupling does not affect the results, since we 
recompute the non-linear parts of eqs. || and [r] after 
each solution of the linearized equations. 

As mentioned before, the constraints are only solved to 
get self-consistent initial data. The constraints are not 
solved during the evolution. 

C. Coordinate Evolution 

In general relativity you have the freedom to choose the 
coordinates however you wish. This freedom manifests 
itself in the choice of original coordinates as well as the 
evolution equations for the lapse a and the shift (3\ It 
is very easy to make a bad choice of coordinates. For 
example, consider a small perturbation in the original 
coordinates x M 

a:" a:" + f. (11) 

A seemingly benign choice of evolution equations for a 
and (3 l can make the small perturbation grow exponen- 
tially, creating coordinate singularities in an otherwise 
non-singular evolution. To prevent any bad behavior in 
£ M , we can use harmonic coordinates. These coordinates 
satisfy 

□a;' 1 = 0. (12) 

That is, the coordinates obey a wave equation (although 
with covariant derivatives). Since eq. |l2| is linear, and 
the original coordinates x M obey eq. [l^, we then get the 
condition on £ p 

= 0. (13) 

Thus, ^ should be wave like, and not exponentially in- 
creasing. 

Imposing eq. [l^ implies evolution equations for the 
lapse and shift. In practice, we use the evolution equation 
for g tt = —a 2 + (3 l f3i instead of the lapse. Then the 
equations are 

^ = ( 7 « a 2 _ pipj) ( _ mt + 2/3 . .) + 2 /?V M (14) 
and 

- (Y J a 2 - ( 7y , fe - 2 lkj>i ) + g tt , k - (15) 



D. Matter Evolution 

Our method is exactly the same as in IlH] except for our 
treatment of low density regions and our time stepping. 
We briefly summarize the method here. 

We represent the matter with three variables — D, t , 
and Si — which roughly correspond to the mass density, 
energy density, and momentum density. We can then 
write the evolution equations in the form 

d (variable) . , 

-^—qI " + ^(flux) 1 = (source), 

where the only spatial derivatives of the matter variables 
occur in the di(flux) 1 term. Writing the equations in this 
form allows us to apply the cornucopia of shock capturing 
methods developed over the years for ordinary hydrody- 
namics. We use a Roe scheme with a standard min-mod 
piecewise-linear reconstruction algorithm |24[ . 

The fluxes and sources are written in terms of prim- 
itive variables {p,e,Vi} which are not simply related to 
the evolved variables {D,r,Si}. Thus we have to do a 
costly root finding at each point at each time level to find 
{p, e, Vi} from {D,r, Si}. This step can go awry in low 
density regions, because numerical errors can conspire 
to evolve {D,t, Si} into something that has no physi- 
cal {p, e, Vi}. To combat this, we use a variation of the 
method of |15|] . We create a fake "atmosphere" in the 
empty space around the star by setting a minimum value 
of D of about 1CP 9 of the initial central density of the 
star. Whenever D is evolved to a value lower than 10 -9 , 
we set D to 10~ 9 . Also, if D is lower than 10~ 5 , we set r 
and Si to zero. Furthermore, if we ever get a transform 
from {D,t, Si} to {p,e,Vi} that gives unphysical values 
(e.g., negative energies), we start over, replacing the defi- 
nition of t with the condition for adiabatic flow P = kp r . 
We then compute {p, e, v{} and recompute r. This simple 
prescription allows long term evolutions of neutron stars 
and their surroundings. 

Apart from the low density treatment, the only other 
difference between us and |T|] is the time stepping. We 
implement a Strang split of the hydrodynamics [p4[ . 
That is, we evolve for a half time step as if the evolution 
equations are only the flux terms of the hydrodynamics. 
Then we evolve for a full time time step as if the evolu- 
tion equations are only the hydrodynamic source terms 
and all the metric and coordinate terms. Then we again 
evolve for a half time step as if the evolution equations 
are only the flux terms. 

When we evolve the fluxes for the first half step, we 
evolve first in one randomly chosen direction d\, then in 
another randomly chosen direction c?2, and then in the 
third direction 0I3. For the second half-step, we do it in 
the reverse order dz-d2-di. Thus, for each step we choose 
from one of six random orders: xyz-zyx, xzy-yzx, yxz- 
zxy, yzx-xzy, zxy-yxz, and zyx-xyz. 

We use second-order Runge-Kutta stepping for the 
matter, metric, and coordinate equations. That is, when 
evolving the equation 
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from t to t + At, we first take an intermediate step 

^intermediate — Qt ~l~ ^ ^\Qt)i 

and then we take the full step 

Qt+At = Qt + At ^(^intermediate)- 

In conjunction with the Strang split, this makes the whole 
evolution second-order accurate. We use a Courant factor 
At /Ax = 0.25/c. This scheme is different from [[TH). 
They evolve the fluxes in the x, y, and z direction all at 
once, and their coupling in time between hydrodynamics 
and spacetime is very different. 



E. Boundary Conditions 

On the outer boundaries, we use the interpolated Som- 
merfeld condition of |l7j for the metric and coordinate 
terms. Essentially, we assume that all of the metric and 
coordinate variables behave like outgoing, radial waves 



Q(t,r) 



G(at - (det7)*r) 



where Q = - 5 l j,K l j,g u + 1, ft}. Thus, we compute 
the value at the boundary by following the characteristic 
back to the previous time step and linearly interpolat- 
ing the corresponding variable to that point. We tried 
just freezing the boundaries (i.e., 7oid = 7ncw), but that 
seemed to reflect and amplify waves, leading to an insta- 
bility. 

For the matter terms, we use "flat" boundaries. That 
is, values at points on the boundary are set equal to those 
at the points just inside the boundary. Although this con- 
dition seems crude, it works rather well at propagating 
blobs of matter off of the computational grid pq] . 



F. Numerical Diffusion 

Evolutions using the methods described up till now 
have failed catastrophically because of short wavelength 
numerical instabilities that grow without bound. To con- 
trol these instabilities, we add some numerical diffusion. 
We do this by adding a V 4 term to the right hand sides of 
eqs. ||, [| 0j and [l5|. This technique, introduced in [ p6[ , 
has the effect of spreading out short wavelength, poorly 
resolved features, while leaving long wavelength, well re- 
solved features alone. This is exactly what is required — 
we don't need an accurate evolution of short wavelength 
features since they are poorly represented on a grid with 
finite spacing. Thus, eq. || becomes 

^ = -2aKij + Vifc + Vjft - q (Ax) 3 V 4 7 y , (16) 



where q is a small constant (we use q — 0.09/i?2, where 
i?* is the radius of the neutron star). The (Aa;) 3 fac- 
tor ensures that, as the resolution increases, the mod- 
ified equation quickly converges to the original contin- 
uum equation. We implement the V 4 as V 2 (V 2 ), and 
assume that V 2 = on the boundaries. This is com- 
patible with the \/r falloff, since V 2 (1/r) = on the 
boundaries. The usual way to implement V 2 is as cen- 
tered second derivatives along the coordinate axes (V 2 = 
d 2 /dx 2 + d 2 /dy 2 + d 2 /dz 2 ). Thus, a plot of the finite dif- 
ference stencil in the x-y plane looks like Fig. [l| 

This gives diffusion that acts primarily along the co- 
ordinate directions x, y, and z. Unfortunately, eq. ^ 
has mixed derivatives of 7y (e.g., d 2 jij/dxdy), so points 
couple along diagonal directions. To fix this, we define 
new, diagonal directions u = x + y,v = x — y, w = z, 
and take derivatives along these directions. Then the fi- 
nite difference stencil looks like Fig. ||, which couples 
along the diagonals in the x-y plane. We repeat this for 
the x-z and y-z planes and average the three different 
representations of V 2 . This gives us a stencil that cou- 
ples along the axes and diagonals. Finally, we apply this 
"diagonal diffusion" to the new time level. That is, or- 
dinarily, without the diffusion, the evolution equation is 
implemented as 



7n 



7old + At(RHS), 



where RHS is the right hand side of eq. 
diffusion by changing the update to 



We add 



7n 



= 7old + At{RHS - qAx 3 V 4 [ lold + At (RHS)}}. 



This is equivalent to first order to eq. We imple- 

ment it this way because short wavelength instabilities 
can sometimes change sign at each time step. Applying 
numerical diffusion to the old time step as in eq. [16] can 
end up being always a step out of sync with the instabili- 
ties, adding when it should be subtracting. Implementing 
it in this way doesn't allow the other parts of the equation 
(RHS) to interfere with the diffusion. 

Adding explicit diffusion does slightly alter the equa- 
tions away from the original physical ones, although the 
change gets smaller with higher resolution. In that sense, 
you can just consider it a slight addition to the error, get- 
ting smaller with more resolution. However, the error is 
not random; it moves error in the direction of too much 
diffusion. In an attempt to minimize this error, we tried 
applying this diffusion only to the coordinate evolution 
equations (eqs. [li], [H]), because the coordinate evolution 
does not affect physical quantities. Yet even when we ex- 
plicitly enforced the constraints at each time step, short 
wavelength spikes grew and destroyed the simulation. 

This diffusion is very effective. Ordinarily, second- 
order Runge-Kutta is unconditionally unstable for most 
problems. In this case, instabilities seem to be domi- 
nated by other factors. We have tried alternative step- 
ping alorithms, such as iterative Crank-Nicholson with 
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varying numbers of iterations, and, in every case, the 
evolution was unstable without the diffusion, and stable 
with it. 



G. Computational Issues 

The entire code is written in CH — h and runs on the 
Cornell Theory Center IBM SP2 supercomputer. The 
code has also been ported to a Sun workstation, and is 
freely available from the authors. The runs with 9 3 , 17 3 , 
and 33 3 points used 1 processor, the 65 3 runs used 8 pro- 
cessors and the 129 3 runs used 64 processors. We used 
the Kelp libraries |27j to split up the grid among multi- 
ple processors and simplify the message passing between 
processors. The code scales from 1 processor to 64 pro- 
cessors with about 70% parallel efficiency. All of the runs 
(except for the convergence tests, which were very short) 
took f 2-72 hours, taking f 000-f 500 steps every f 2 hours. 
The hydrodynamics used about 80% of the CPU time, 
the metric and coordinate evolution used about 15%, and 
everything else took up the remaining 5%. 

III. TESTS 

We performed two different types of tests: short-term 
tests to verify the correct implementation of the equa- 
tions, and long-term tests to verify the stability and ac- 
curacy of the method. 

A. Short-Term Tests 

To test the hydrodynamic evolution and its ability to 
handle shocks, we ran shock tubes as in p5| . We set up 
a one dimensional shock tube with p = 10, P = 13.3, 
v = 0, on the left and p = 1, P = 0.66 • 10~ 6 , v = on 
the right. This problem has an analytic solution [ ^8|j29| 
which we can compare to our evolution. Figure |^ shows 
a comparison of the normalized values of the pressure, 
density, and velocity, and we can see that our method 
reproduces the analytic answer reasonably well. 

To test the metric and coordinate evolution, we ran 
convergence tests with the exterior metric of a black hole 
in full harmonic gauge |30| j. To test all of the evolu- 
tion terms and their coupling to each other all at once, 
we ran convergence tests on static and boosted neutron 
stars. We performed our neutron star tests with a T = |, 

K = 5.380 • 10 9 cm 4 g _2 / 3 s~ 2 polytrope with a central 
density /d CC ntrai = 1 4 10 15 gcm~ 3 , resulting in a mass 
M* = 1.35 • 10 33 g, a Schwarzschild radius R = 12.7km, 
and thus a compaction M/R = 0.08. The maximum cen- 
tral density for this equation of state is about /? max = 
3.79 • 10 15 gcm~ 3 , corresponding to a maximum mass of 
1.51 • 10 33 g. Then we converted the variables into the 
harmonic gauge f30| which gives us a harmonic radius 



of i?* = 11.7km. We boosted the stars with a variety 
of velocities [if = (0.3c, 0.2c, 0.1c), (0.9c, 0, 0), (0, 0.5c, 0)] 
using the prescription in jn|. We did convergence tests 
over both a cube entirely within the star and a cube en- 
tirely outside the star. We offset the grid by 10 -5 i?* in 
the x, y, and z directions to avoid having to deal with 
spherical coordinate pathologies when setting up the ini- 
tial configuration. 

When doing convergence tests, we first took one step 
on a 9 3 grid, two steps on a 17 3 grid, and four steps on 
a 33 3 grid. All three grids spanned the same space, so 
each step doubled the resolution. Then we compared the 
solution on the three grids against the analytic solution. 
Because our differencing and time stepping are both sec- 
ond order, the error should decrease as the square of the 
number of points (i.e., a factor of four with each step in 
resolution) . 

However, this is not a perfect argument because there 
are higher order errors contributing to the overall error. 
Fortunately, these errors decrease even faster than the 
second order error. Therefore, if we plot the 9 3 error, the 
17 3 error multiplied by 4, and the 33 3 error multiplied 
by 16, then we should see the 17 3 error closer to the 33 3 
error than the 9 3 error. 

Unfortunately, even this procedure is not foolproof. 
Sometimes the error can go through zero, leading to cases 
where the 17 3 error is closer to the 9 3 error than the 33 3 
error in isolated regions. Figure [| shows a plot where 
this happens. Sometimes the higher order error is still 
large, leading to cases where the high resolution error 
decreases so quickly that the medium resolution error is 
no longer closer to it. As an added difficulty, our outer 
boundaries aren't very good, so we couldn't use the error 
near it. In practice, we had to throw away the two outer- 
most points, leaving only 5 3 = 125 points to compare. In 
addition, our solution of the initial data for the neutron 
star is imperfect, so the center point was often bad (For 
example, the only bad point for j xy is at the center). 

Finally, to get proper convergence, we had to turn off 
the numerical diffusion terms, because they add unphysi- 
cal, although small, diffusion, and also the min-mod slope 
limiter in the hydrodynamics, because it enforces first or- 
der convergence near maxima. Keeping these techniques 
allows longer evolutions while sacrificing strict second- 
order convergence. 

Even with all of these caveats, checking convergence 
was a powerful tool for finding bugs. An error in how we 
implement an equation prevents the entire grid from con- 
verging. Thus, a plot like Fig. [|, which shows the worst 
converging variable for that particular set of convergence 
tests, was not worrisome. We can trace its defects back 
to a combination of what is wrong in Fig. ^ and prob- 
lems at the center. When we can not do that, the cause 
is a bug. We found many bugs this way. For example, 
we accidentally implemented the —2aKnK l j term in eq. 
|3| as —2aKuK 3 l7 and the indexing error prevented con- 
vergence. This sensitivity gives us a reasonable degree of 
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confidence that we implemented the correct equations, 
and implemented them without error. 



B. Long-Term Tests 

To test the ability of our code to handle long neutron 
star evolutions, we ran the code with a single, stationary 
star. The idea is to see how well it holds a static star 
in equilibrium. In this section and later we adopt the 
units R* = 11.7km = 1, M* = 1.35 • 10 33 gm = 1, and 
c = 1. This implies G = 0.08, 0.039ms = 1, and the 

— 1/2 

hydrodynamic timescale p ccntra i = 4.49. We compute 
the norm of a variable by averaging its square over the 
grid and then taking a square root. 

The first set of tests evolved just the metric and co- 
ordinate terms, keeping the matter terms fixed to their 
analytic values. Figure ^| plots the error in the energy 
constraint (eq. ||) vs. time for three cases. It quickly 
drops down and levels off as the simulation settles into 
a stable configuration. Ideally, we would like that con- 
figuration to be the analytic one, but finite resolution, 
numerical diffusion, and the imperfect outer boundary 
condition will change it. Figure plots a cross section 
of the error for the three cases in the final state. The 
errors for the low resolution, close and far boundaries 
cases are nearly identical. The only visible difference is 
near x = 2, where the close boundary has a little more 
error because of the boundary. It is also clear that bet- 
ter resolution helps a lot, drastically decreasing the error. 
This is probably because the biggest contributor to the 
error here is the numerical diffusion terms. Their effect 
decreases as the cube of the grid spacing, which is much 
faster than the other errors, which decrease as the square 
of the grid spacing. The reason that the constraint norm 
for the far boundary case is so much smaller than the 
close boundary case is that the error is concentrated in 
the center near the star. The far boundary case just has 
more points over which to average the error. This is en- 
couraging, because it means our simulation will improve 
as we increase resolution (not always a given), and our 
boundary is probably adequate, even relatively close in. 
Figure || plots a cross section of the relative error in 
det "fij . Even for the worst case, the errors are only 



about 1%. The effect of the boundaries is much stronger 
than in Fig. ^. The simulation does not stray far from the 
analytic solution, giving us hope that numerical diffusion 
and the boundary conditions do not significantly affect 
the results. 

The second set of tests evolved just the hydrodynamics, 
keeping the metric and coordinate terms fixed to their 
analytic initial values. Figures || and |l0| show a cross 
section of Tff versus time for two resolutions. The star 
rapidly diffuses out until it starts to interact with the 
boundaries. Remember that the hydrodynamic terms do 
not have any numerical diffusion. This diffusion doesn't 
actually cause a problem until late times, when "fingers" 



appear, stretching from the star to the boundary (Fig. 
|l l| ) . These fingers quickly grow until they dominate the 
simulation. The simulation still looks well behaved, if 
inaccurate, for a long time before then (Fig. [T2|) . 

The third set of tests evolved everything: hydrody- 
namics, metric, and coordinates. The results are largely 
the same as the hydrodynamics-only results (Figs, [l^, 
|l4| , and |l5|). We also ran a test with far boundaries, 
and, as expected, we were able to evolve for significantly 
longer times. Figure |l6| shows a plot of the central den- 
sity T£ vs. time for the three cases. Eventually, the 
star again developed "fingers", at which time the code 
became hopelessly inaccurate. 



IV. COALESCENCE 

To show that the method will work for its intended 
purpose, simulating coalescing stars, we simulated a coa- 
lescence with easily computed, although not astrophysi- 
cal, initial data. We took the equilibrium solution for the 
neutron star, replicated it next to itself, and then boosted 
the stars with v = 0.15c in opposite directions as shown 



in Fig. 17 (the Kepler frequency for two point particles 
is v — 0.19c). To get self-consistent gravitational initial 
data, we solved the constraints with this matter source. 

This is not astrophysically interesting initial data by 
any means. The equation of state, and thus the size, 
mass, and internal structure of the stars, are all wrong. 
Even if they were right, setting up the variables by plac- 
ing boosted solutions next to each other is definitely 
wrong. There is likely to be a large amount of initial 
wave content, which can seriously affect the dynamics. 
We discovered this to our chagrin when we tried starting 
the stars a little farther apart. The initial wave content 
pushed the stars away from each other! 

In addition, the resolution of these te sts is really no 
better than the long term tests in section [II B, which, as 
we saw, were not very accurate. Even so, it serves as an 
important validation of the method. 

To measure the gravitational waves, we adopt a scheme 
similar to that used in |3lj . We can define the transverse 
traceless (TT) gauge wave amplitudes h + and h x by pro- 
jecting out components of the metric. Along the z axis, 
this becomes 



h% — Ixy 



Ivv) 



(17) 
(18) 



A more exact method would be to use gauge-invariant 
Moncrief variables J32j and integrate over spherical har- 
monics p3| |. Considering the accuracy of the underlying 
simulation, this would be overkill. 

Figure [l^ shows the gauge invariant trace T£ vs. time 
for a low resolution run with the boundaries very far out 
at (—8,8). The stars complete about two orbits before 
completely merging. 
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Figures |l9| and gO| show the amplitude of h + and h x 
for three resolutions. As we increase the resolution, the 
waveforms seem to converge, revealing finer and finer de- 
tails. The atmosphere was treated a little differently in 
these runs, with the floor value of D set at 10 _4 D ce ntrai 
instead of lCP 9 -D ccn trai- This causes the evolution to go 
bad earlier, around t = 60. The earlier waveforms are 
not significantly affected, though. 

If we run some simulations with the boundaries farther 
out, the results are not as encouraging. Plotting two res- 
olutions of far boundary runs with the close boundary, 
high resolution run (Figs. |2l], we see that the initial 
amplitude and arrival time of the h x wave is different. It 
arrives later, whereas if it were arriving from the z ~ 
plane, where the stars are, it should arrive at an earlier 
time. This can not be explained from differences in how 
the waves propagate from r — 4i?* to r — 8i?* . Figure |24| 
shows wave amplitudes for the low resolution, far bound- 
ary run extracted at different radii. The h x amplitudes 
match each other very well. 

Looking at the h + wave, we see new structure around 
t = 10—15 (twin spikes) that is absent in the close bound- 
ary runs. This is probably just propagation effects, for if 
we plot the wave amplitudes extracted at different radii 
(Fig. 23), we see that they match the structure much 
more at z = 4. 

These differences in h x probably arise from differences 
in the initial data. Remember that we solve the con- 
straints to get the initial data. We assume that the outer 
boundaries are correct, and change the interior metric 
variables to fit. Thus, the different simulation domains 
will have varying amounts of initial wave content. This 
also suggests, unfortunately, that the big, initial peaks in 
h+ and h x may not be real. 

The stress energy tensor for weak waves in the TT 
gauge is @ 



where the (...) symbol means an average over several 
wavelengths. The total energy that passes through a thin 
shell at radius r in time At is 4:TrT 00 r 2 cAt. To estimate 
hijfi, note that h+ and h x go from ~ 10~ 2 to over a 



time of ~ 10-R*. We extract the waves at 4i?* 
luminosity is 



so the 



L ~ 4tt 



1 (h 2 + + h 2 x ) 

2 (10i?*) 2 



2 • 10" 



If we compute an estimate of the gravitational wave 
luminosity from the quadrupole formula for point masses 
in circular orbits |55L Section 36.6], 

L=^^!~3-10- 
5 a 5 

where [i = mimi/M, M = mi + rri2, a = separation = 
2R* . The large differences in magnitude suggest that the 



waves may be generated by something other the stars, 
such as the interaction between the right and left sides of 
the grid (Fig. [It]). The spacetimes are initially boosted 
in opposite directions, so there is a discontinuity at x = 
0. This discontinuity is smoothed somewhat, but not 
removed, when we solve the constraints. 



V. FUTURE DIRECTIONS 

There are a few things that must be done before the 
code will produce interesting results. First, the code 
needs realistic initial data along with more realistic equa- 
tions of state. A number of groups |56]-[58| have managed 
to construct reasonable initial data. Getting these results 
and converting them into a usable format for the dynam- 
ical code is, in principle, a straightforward task. 

Second, the code needs better resolution. Simply in- 
creasing the number of points is impractical, since the 
code already taxes the capabilities of current supercom- 
puters. Instead, the solution is probably adaptive mesh 
refinement (AMR) 39 1 . Massively parallel machines will 
still be required, but at least then the problem is possi- 
ble. Unfortunately, AMR on massively parallel machines 
is still a new endeavor, so libraries for making it straight- 
forward to change a normal, parallel code into parallel 
AMR code are still immature. AMR will also allow us 
to move the boundaries out much farther, relieving some 
headaches there. 

These steps will allow the code to be applied to the 
most pressing reason for these simulations: providing ac- 
curate theoretical waveforms for the new generation of 
gravitational wave detectors coming on line. It will also 
give reasonably accurate estimates of the amount of ma- 
terial ejected as possible r-process elements. To under- 
stand 7-ray bursts, on the other hand, will require this 
and more. Implementing magneto-hydrodynamics will 
allow us to evaluate, for example, the likelihood of the 
DRACO model jl(|. Adding neutrino generation and 
transport will account for an important energy source. 
Also, neutrino effects may, as in supernovae, affect the 
dynamics El. 
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FIG. 1. Finite difference stencil along axes. The circles 
denote which points V 2 samples in the x-y plane. The points 
with X's do not affect the value of V 2 . 



FIG. 2. Finite difference stencil along diagonals. 



FIG. 3. Normalized values of the density p, pressure P, 
and velocity v for a one-dimensional shock tube at t — 0.5 
using 256 grid points. 



FIG. 4. Plots of the sign of the difference 

|4 • | IT 3 error] - |9 3 error] | — 1 16 ■ 1 33 3 error] - 4 ■ |l7 3 error] | 
for the variable D for a neutron star with the boost 
if — (0, 0.5c, 0). The boundaries are at ±0.257?*, but the 
outermost 2 points are discarded, so this grid is only 5 points 
on a side. When the difference is positive, the normalized 
17 3 error is closer to the normalized 33 3 error than the 9 3 
error, implying correct convergence. Thus, negative numbers 
(the dark areas) indicate bad convergence. Note that the plot 
only goes negative in limited regions. The 3rd order error 
goes through zero there, so the normalized errors cross each 
other. 



FIG. 5. Plots as in Fig. ^ but for the variable K yz . For 
this set of convergence tests, this is the variable that looks 
the worst. 







FIG. 6. The norm of the error in eq. g vs. time when the 
matter terms are kept frozen for three cases: a low resolution 
run, a high resolution run, and another low resolution run 
with the boundaries twice as far out. 

FIG. 7. The error in eq. ^ for the relaxed state on Fig. ^[ 
The star is at the origin, and the star is spherically symmetric. 
This plot is taken from the x > 0, y — z — line. 



FIG. 8. A cross section of the relative error in */ det 7^ 
for the relaxed states of Fig. [| 

FIG. 9. A cross section of Tj? for a low resolution run (33 3 
points) when only the hydrodynamics are evolved. The star 
quickly spreads out, although the scale is very small. 



FIG. 18. Tff for coalescing stars. The simulation spans 
(—8,8), but we only plot (—2,2) to show the details of the 
merger. 

FIG. 19. h+ (eq. 0) for three different resolutions ex- 
tracted at the boundary z — 4. 

FIG. 20. h x (eq. [if) for three different resolutions ex- 
tracted at the boundary z — 4. 

FIG. 21. h+ (eq. 0) 

for the high resolution run of Fig. 
[Hj , compared to low and medium resolution runs with the 
boundaries twice as far out. The farther boundary wave am- 
plitudes are extracted at z = 8 and scaled to account for the 
1/r falloff. The structure differs markedly, especially around 
t = 10 - 15. 



FIG. 10. A cross section of for a high resolution run 
(65 3 points) when only the hydrodynamics are evolved. The 
star spreads out much more slowly than the low resolution 
run of Fig. |{)| 

FIG. 11. A cross section of in the equatorial plane 
taken at time t = 68.5. A finger has developed along the 
x = line. 

FIG. 12. A cross section of T£ in the equatorial plane 
taken at time t — 50. There is no evidence of a finger. 



FIG. 22. h x (eq. Q for the high resolution run of Fig. 
, compared to low and medium resolution runs with the 
boundaries twice as far out. The farther boundary wave am- 
plitudes are extracted at z — 8 and scaled to account for 
the 1/r falloff. The initial peak is significantly stronger and 
arrives later. 

FIG. 23. h+ for a low resolution run (65 3 points) with 
the boundaries at (—8,8). The amplitudes are extracted at 
z — 4 and z = 8. They do not match up, suggesting that the 
differences in Fig. El] may be due to propagation effects. 



FIG. 13. A cross section of for a low resolution run 
(33 3 points) when everything is evolved. Note that it is very 
similar to Fig. 

FIG. 14. A cross section of for a high resolution run 
(65 3 points) when everything is evolved. Once again, very 
similar to the hydrodynamics-only result (Fig. 



FIG. 24. h x for a low resolution run (65 3 points) with 
the boundaries at (—8,8). The amplitudes are extracted at 
z — 4 and 2 = 8. They match up reasonably well, especially 
the large initial peak. 



FIG. 15. A cross section of T£ for a low resolution run, 
but with the boundaries twice as far as in Fig. Note that 
it takes much longer to reach the boundaries, resulting in a 
much longer evolution. 



FIG. 16. The central density of the star for three different 
runs. Note that the two low resolution cases track each other 
very well until the close boundary case diverges dramatically. 



FIG. 17. Overhead view of the initial state of the two neu- 
tron stars 
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